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Experiments on fracture surface morphologies offer increasing amounts of data that can be ana- 
lyzed using methods of statistical physics. One finds scaling exponents associated with correlation 
and structure functions, indicating a rich phenomenology of anomalous scaling. We argue that tra- 
ditional models of fracture fail to reproduce this rich phenomenology and new ideas and concepts 
are called for. We present some recent models that introduce the effects of deviations from ho- 
mogeneous linear elasticity theory on the morphology of fracture surfaces, succeeding to reproduce 
the multiscaling phenomenology at least in l-fl dimensions. For surfaces in 2-1-1 dimensions we 
introduce novel methods of analysis based on projecting the data on the irreducible representations 
of the SO(2) symmetry group. It appears that this approach organizes effectively the rich scaling 
properties. We end up with the proposition of new experiments in which the rotational symmetry 
is not broken, such that the scaling properties should be particularly simple. 



It is a privilege to dedicate this paper to Pierre Hohenherg and Jim Langer who contributed, 
respectively, decisive ideas to scaling concepts in statistical physics and to the study of fracture. It 
is our hope that the marriage of these two issues in the present paper may give them some degree of 
pleasure. 



I. INTRODUCTION 

The failure of materials is a challenging interdisci- 
plinary problem of both technological and fundamental 
interests. From the technological point of view, the un- 
derstanding of the failure mechanisms of materials under 
various external conditions may improve dramatically the 
integrity of structures in a wide range of applications. 
From the theoretical point of view, the understanding 
of the way materials fail entails the development of new 
mathematical methodologies and necessitates the intro- 
duction of new concepts in non-linear and solid state 
physics. The modern development of the field as a sci- 
entific discipline initiated with the pioneering work of 
Griffith [1] who identified the importance of defects in 
determining the strength of materials. These defects act 
as stress concentrators in the sense that the typical stress 
near the defect can be much higher than the applied 
stress. In that way the strength of materials is highly re- 
duced, explaining the long lived conflict between theoret- 
ical strength estimations and experimental observations. 
In order to understand the way defects affect the failure 
of materials we first introduce a simple scaling argument 
in the framework of equilibrium thermodynamics. 

In Fig. 1 we consider a material under the application 
of a uniform external stress tioo at its far edges. In addi- 
tion, the material contains a single defect that is assumed 
here to be a crack of length L that is cutting through the 
material. A crack is a region whose boundaries cannot 
support stress. In the absence of a crack the material 
is uniformly stressed with an associated potential energy 
density £p 

fp-#, (1) 




FIG. 1: A material under the application of a uniform exter- 
nal stress (Too at its far edges in the presence of a crack of 
length L that is cutting through the sample. The shaded re- 
gion represents the typical area in which the potential energy 
density is changed relative to the uniform stress state. 

where E is Young's modulus and it is assumed that the 
material is linear elastic, resulting in a second-order en- 
ergy density. The presence of a crack of length L releases 
the stresses in an area of the order of ^ (shaded area 
in Fig. 1), resulting in a reduction AUp in the potential 
energy per unit material width 

AC/p = -cp^ , (2) 
where Cp is a dimensionless factor. On the other hand. 
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the very existence of the crack is associated with free 
material surface. Assuming that the energy cost per unit 
area of free surfaces is F, tlic generation of a crack of 
length L results in an increase in the surface energy per 
unit width by an amount Af/^, 



A[/s = TL 



The total energy change per unit width is 



At/ = -c„ 



E 



+ TL . 



(3) 



(4) 



This equation tells us that for small values of L the for- 
mation of a crack is costly (AC/ > 0) whereas longer 
cracks are energetically favorable (AC/ < 0). Actually, 
once a critical length is achieved the crack tends to in- 
crease indefinitely until the material completely fails. 
This non-equilibrium catastrophic crack propagation is 
at the essence of the failure of material. The two differ- 
ent regimes described above are separated by a Griffith 
critical length 
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which is shown to be a combination of material properties 
[E and F) and external loading conditions {(Jqc)- 

The thermodynamic argument predicts that crack 
propagation should always be catastrophic once initi- 
ated. In fact, although fast crack propagation is com- 
mon, there are many situations in which the crack evolves 
quasi-statically. The point to stress is that whatever is 
the mode of crack growth, the argument exemplifies the 
multi-scale nature of the phenomenon; the potential en- 
ergy released from the large scales dissipates in a very 
localized region near the crack tip where new crack sur- 
faces are generated. 

The technical discussion of crack propagation is clas- 
sically done in the context of "linear elasticity fracture 
mechanics" . The state of deformation of the material is 
described by the displacement field u(r). This field con- 
sists of three translational degrees of freedom, where the 
local rotational degrees of freedom are neglected [2]. To 
develop a theory that is both translational and rotational 
invariant one defines the strain tensor e^, which is the 
symmetric part of the gradient of u(r) 



{diUj + djUi) 



(6) 



In a material that is homogeneous and isotropic, only 
the invariants of e^j can appear in the expression for the 
strain energy density £. Restricting the analysis to a 
linear theory, i.e. to a quadratic strain energy density, 
we arrive at 
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(2MTr(e2)+ATr2(e)) , 



(7) 



where the Lame coefficients fi and A are material param- 
eters that are related to the more common engineering 
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FIG. 2: The typical symmetry modes of fracture, see text for 
more details. 



constants E (Young's modulus) and v (Poisson's ratio). 
The stress tensor cry is related to the strain tensor via 



dE 
den 



which leads to a linear stress-strain relation 



(Jij = 2ij,€ij + XSijCkk ■ 



(8) 



(9) 



The stress field is just the local force per unit area. 

Therefore, the Newton's equations of motion for a unit 
volume of mass density p are given by [3] 



= P- 



(10) 



Substituting Eq. (9) in the last set of equations we arrive 
at the Lame equation 



{X + li)V{V ■ u) + iJ.V^u = p 



9^1 



(11) 



Up to now the field equations for the medium were 
constructed disregarding the effects induced by the pres- 
ence of a crack. At the level of linear elasticity fracture 
mechanics, the crack is introduced as boundary condi- 
tions. As was mentioned before, a crack is a region whose 
boundaries are free surfaces. Denote the unit normal to 
the crack surface at any point by fi; the boundary con- 
ditions on the crack surface are 



aij Uj 



. 



(12) 



These boundary conditions introduce non-linearities into 
the problem even if the field equations themselves are lin- 
ear. This is a major source of mathcmatic;al difficulties. 

It is conventional to decompose the stress field under 
general loading conditions to three symmetry modes with 
respect to the fracture plane. These are illustrated in Fig. 
2. In mode I the crack faces are displaced symmetrically 
in the normal direction relative to the fracture xz plane, 
by tension. In mode II, the crack faces arc displaced anti- 
symmetrically relative to the fracture xz plane, in the x 
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direction, by shear. In these two modes of fracture the 
deformation is in-plane (xy) . In mode III, the crack faces 
are displaced anti-symmetrically relative to the fracture 
plane, in the z direction, by shear. This is an out-of-plane 
fracture mode. 

The dissipation involved in the crack growth, quan- 
tified by the phenomenological material function F, is 
assumed to be highly localized near the crack front. 
Therefore, one is interested in the near crack front fields. 
Within linear elasticity theory these are actually singu- 
lar. To see this consider first a quasi-static infinitesimal 
extension SL of the crack. Denoting the stress field near 
the tip by a(r), the energy (per unit width) released from 
the linear elastic medium dU is 



6U 



SL 9 



E 



dr 



(13) 



This amount of energy is invested in creating new crack 
surfaces whose energy cost (per unit width) is TSL. 
Therefore, we must have 



£7(r) 
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(14) 



The inverse square-root singularity seen in Eq. (14) ex- 
ists also in the fully dynamic case. In fact, by asymptotic 
expansion of the local stress tensor field near the crack 
front, it can be shown that the leading term is given by 
a sum of three contributions corresponding to the three 
modes of fracture [4] 

aij{r,e,t) = 



{9,v) ,.^ii{0,v) (e,w) 



(15) 

where (r, 6) are local polar coordinates system, t is time, 
v(t) is the local instantaneous front velocity, S^*^(^,t;) 
are universal and Ki^^-j are the stress intensity factors. 
The stress intensity factors are non-universal function- 
als of the loading conditions, sample geometry and crack 
history. The predicted singular behavior is, of course, 
not physical and there must exist mechanisms to cut off 
this apparent singularity. Nevertheless, there are many 
physical situations for which the size of the region where 
linear elasticity breaks down is small compared to other 
relevant lengths. Therefore, the stress intensity factors 
are very important physical quantities and the singular- 
ity may be retained in many models. This singularity is 
a source of mathematical difficulties and physical riddles. 
The stress concentration quantified by the stress intensity 
factors shows that the material near a crack front expe- 
riences extreme conditions; the response to these condi- 
tions is far from being well-understood. 

Up to now we have considered the case of a sin- 
gle straight crack in an otherwise linear homogeneous 
medium. A more realistic description of materials in- 
cludes other sources of heterogeneity which exist in every 



material at some level. The effect of distributed soinces 
of "disorder" calls for a statistical treatment. In devel- 
oping such an approach one should understand the in- 
terplay between various fields and the material disorder, 
resulting in a complex spatio-temporal behavior. This 
complexity was studied extensively in the last decade, 
accumulating a wealth of relatively high precision exper- 
imental data and offering a real challenge for the theoret- 
ical physicist [5] . In this paper we focus on one important 
aspect of the problem: the statistical physics of the mor- 
phology of quasi-static fracture surfaces. This subject 
has attracted a lot of interest recently and became a very 
active field of research [6] . The pioneering experimental 
work described in Ref. [7] drew attention to the fact that 
fracture surfaces are rough graphs in 2+1 (1+1) dimen- 
sions when the broken sample is three dimensional (two 
dimensional) and therefore might have anisotropic scal- 
ing properties. Suppose that the surface is described by 
its height h{x, y) at position (x, y) in a smooth reference 
plane. Consider now the probability density P{Ah, i) 
that the height difference between two points in the refer- 
ence plane separated by a length I is within d{Ah) of the 
height difference Aft. The anisotropic scaling properties 
of the surface manifest themselves through the following 
invariance under affine scale transformation 



X"P{X"Ah,Xi) = P{Ah,e) 



(16) 



where H is the roughness exponent. This is the starting 
point for the discussion of self affine properties of fracture 
surfaces. Later we will see that this definition is too lim- 
ited and cannot capture the rich complexity exhibited by 
fracture surface morphology. The first part of this paper 
focuses on two dimensional fracture where the generated 
rupture lines are 1+1 dimensional graphs. In Sec. II 
we show that the statistical properties of rupture lines 
cannot be fully characterized by the scaling invariance 
of Eq. (16). In that case, the more complex structure of 
the probability distribution function leads to multiscaling 
in contrast with monoscaling implied by Eq. (16). We 
emphasize the properties of 1+1 dimensional disordered 
fracture needed to be explained by a proper theory. Sec. 
Ill offers a short critical review of existing theoretical ap- 
proaches to the problem. We discuss the limitations of 
these approaches and explain why, in our opinion, they 
do not provide a satisfactory description of the underly- 
ing physics. In Sec. IV we describe a new theoretical 
model for the growth of a crack in a two dimensional 
medium in the presence of material disorder. We elabo- 
rate on the mathematical foundations of the model based 
on a recent development in which the method of iterated 
conformal maps was applied to the problem of elasticity 
in the presence of irregular crack geometries. We sum- 
marize the results of the model and show that they meet 
the basic requirements of Sec. II. 

The second part of the paper discusses three dimen- 
sional fracture where the generated surfaces are 2+1 di- 
mensional graphs. In Sec. V we show how the scal- 
ing properties of fracture surfaces can be rationalized 
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by decomposing the height-height structure function into 
the irreducible representations of the SO{2) symmetry 
group. This method offers a new way of understanding 
the anisotropic properties of fracture surfaces in the plane 
of fracture. We propose new experiments in which the 
rotational symmetry is not broken such that the scaling 
properties should be particularly simple. Sec. VI offers 
a summary and outlines for future research directions. 

II. MULTISCALING IN 1+1 DIMENSIONAL 
FRACTURE 

In 1+1 dimensions one denotes the graph as h{x), 
where h is the height of the surface at point x relative 
to a smooth reference line and considers the n*'' order 
structure function Sn{^), 

Sn{e) = {\h{x + e) - h{x)n , (17) 

where angular brackets denote an average over all x. 
These quantities are invariant under affine scale trans- 
formations if they are homogeneous functions of their 
arguments 

SniXe) ^ X'^'"' Sn{£) . (18) 

If C*-"-* is linearly related to n the scaling properties of 
the graph are called "normal" and the graph is statisti- 
cally self-afhne. On the other hand, if C*^"^ is a non-linear 
function of n, the structure functions are multiscaling (or 
anomalous) and the graph under consideration is called 
"multiafhne" . To our best knowledge, all the studies re- 
garding fracture in 1+1 dimensions treated the rupture 
lines as normal graphs, implying that C'"' = where 
H is the roughness (Hurst) exponent. Note that H > 1/2 
indicates the existence of positive correlations in the pro- 
cess generating the surface which implies that an upward 
incremental deviation (relative to the smooth reference 
line) of the rupture line is more likely to be followed 
by an upward deviation than an downward one and vice 
versa [8]. This feature shows that the higher H is, the 
smoother is the surface. Experimental analyses of rup- 
ture lines in quasi two-dimensional materials yielded a 
roughness exponent H whose numerical value was close 
to 0.67 [9-12], suggesting some universality of the sur- 
face generating process. The proximity of this numeri- 
cal value to the exact ratio H — 2/3, characterizing the 
roughness of directed polymers in random media [13], has 
lead some authors to suggest that the two problems are 
in the same universality class [10, 13-15]. In this view, 
fracture is considered as a global minimization problem. 

In a recent work [16] we have shown that the phe- 
nomenology is much richer and fracture lines are multi- 
scaling. An example that provides us with information of 
sufficient accuracy to establish the multiscaling charac- 
teristics is rupture lines in paper. The data acquired by 
Santucci et al. [17] was obtained in experiments where 
centrally notched sheets of fax paper were fractured by 



standard tensile testing machine. Four resulting crack 
profiles h{x) were digitized. Each digitization contained 
a few thousand points, where care was taken to insure 
that the smallest separation between points in x is larger 
than the typical fiber width; this is important to avoid 
the artificial introduction of overhangs that destroy the 
graph property. 
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FIG. 3: A log-log plot of S3{£) as a function of £. The linear 
fit corresponds to a typical scaling range of about 1.5 orders 
of magnitude, with a slope of ~ 1.65. 

Denoting Ah{£) = h{x + £) — h{x) we analyzed the 
data by boxing £ in logarithmic boxes, accumulating the 
data between 10° and 10°'^^ (the smallest box) and be- 
tween 10^-^^ and 10^-^ (the largest box). The m*'* box 
was considered as representing data for I = iO"'^*'-25. 
On the basis of this boxing we constructed the proba- 
bility distribution function (pdf) P(A/i(^)), which is just 
the probability distribution function defined in Eq. (16), 
by combining data from all the four samples. Samples 
that exhibit marked trends (probably due to the finite 
size of the sample), were detrended by subtracting the 
mean from each distribution. The computed pdf's were 
then used to compute the moments (17), and these in 
turn, once presented as log-log plots, yield the scaling 
exponents C*^"). Such a typical log-log plot is shown in 
Fig. 3, exhibiting a typical scaling range of more than 1.5 
orders of magnitude. The resulting values of the scaling 
exponents C''"'' are shown in Fig. 4. 

As a function of n these numbers can be fitted to the 
quadratic function C^"^ = nH-n^X with H = 0.64 + 0.03 
and A = 0.026 + 0.004 (a hnear plot nC^^^ is added for 
reference). The error bars quoted here reflect both the 
variance between different samples and the fit quality. 
The exponents were computed for n < 8 since for higher 
moments the discrete version of the integral 

j \Ahie)\''p{Ah{i))dAh{e) (19) 

did not converge. On the other hand, the convergence of 
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FIG. 4: The spectrum (^'"^ as a function of the moment order 
n for rupture hues in paper. The function is fitted to the form 
(^^"•' = nH — n^X and the parameters H and A are given. The 
errors in the estimation of these parameters reflect both the 
variance between different samples and the fit quality. The 
linear plot n^'^' is added to stress the non- linear nature of 
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FIG. 6: The natural logarithm of P{Ah{e)y" as a function 
of Ah{£)/£" for H = 0.64. The legend gives the scale e for 
each plot. 
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the 8 order moment is demonstrated in Fig. 5. 
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FIG. 5: An example of the convergence of the integral in Eq. 
(19) for e = 10^-^^ and n = 8. 

The point to stress is that the scaling exponents C''"'' 
depend non-linearly on n; for the range of n values for 
which the moments converge, the exponents can be fit- 
ted to a quadratic function. It is well known from other 
areas of nonlinear physics, and turbulence in particular 
[18, 19, 55], that such phenomena of multiscahng are as- 
sociated with pdf's on different scales i that cannot be 
collapsed by simple rescaling. In other words, in the ab- 
sence of multiscaling, there exists a single scaling expo- 
nent H with the help of which one can rescale the pdf's 



p{Ah{e)) 



(20) 



where /(•) is a scaling function. In our case such rescal- 
ing does not result in data collapse. In Fig. 6 the nat- 
ural logarithm of P{Ah{£))£^ is plotted as a function of 
Ah{e)/e" for H = 0.64. Indeed, the data does not col- 
lapse onto a single curve. The fatter tails of the probabil- 
ity distribution functions at smaller scales are typical to 
multiscaling situations, where the statistics of rare events 
plays an important role. 

The multiscaling property of 1+1 dimensional fracture 
surfaces has an important implication; Since the prob- 
lem of directed polymers in random media is monoscal- 
ing [21], we conclude that the two problems are not in 
the same universality class. In fact none of the known 
kinetic growth model of statistical physics [13] is consis- 
tent with the multiscaling spectrum discovered in 1+1 
dimensional fracture. This might suggest that 1+1 di- 
mensional fracture defines a distinct universality class in 
statistical physics, a proposition that must be checked 
against a much broader experimental data set. To con- 
clude, the presence of multiscaling offers a stringent test 
for any theoretical model of 1+1 dimensional in-plane 
fracture. 



III. TRADITIONAL THEORETICAL 
APPROACHES 

In the last fifteen years or so several statistical physics 
approaches were applied in this field. In this section we 
offer a short critical review of these approaches. 
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A. Lattice Models 

The most studied model aiming at understanding the 
origin of self-affinity of fracture surfaces morphology, the 
value of the roughness exponent H and its degree of uni- 
versality is the quasi-static random fuse model [22] . This 
model consists of an array of identically conducting lin- 
ear (Ohmic) resistors placed between two electrical bus 
bars with a fixed voltage difference. A threshold value tj 
taken from an uncorrelated distribution p{t) is assigned 
to each resistor. At each step of the numerical simula- 
tion the Kirchhoff equations are solved and the current 
ij passing through each resistor is calculated. The sys- 
tem is advanced in time according to the extremal rule 
that states that the resistor with the largest ratio be- 
tween current and threshold maxj{ij/tj) is irreversibly 
removed. The system is driven to its final state where 
the network of fuses stops conducting. The behavior of 
this model is controlled by the width of the threshold 
distribution p{t) [23]. For "very strong" disorder, i.e. 
for broad enough threshold distributions, the model be- 
haves like a percolation problem, without any localized 
failure mode. For this regime, the disorder dominates the 
correlations which are induced by current enhancements 
due to removed resistors (the analog of stress concen- 
tration in the elastic problem). The cluster of broken 
bonds (removed resistors) shows no anisotropic scaling 
(formally, the roughness exponent H is one in this case). 
For smaller disorders damage is distributed throughout 
the system, but at the "critical phase" of the dynamics 
of the model a localized failure mode sets in, resulting in 
a self-affine rupture line. At even smaller disorders, in 
the "weak" disorder regime, the system is "critical" from 
the early stages of the evolution in the sense that a large 
crack quickly nucleates and the dynamics is controlled 
by its propagation. The roughness exponent H for the 
regimes where self-affine rupture lines are generated was 
measured to be i? « 0.7 [14, 24]. 

Although the random fuse model is widely regarded "a 
paradigm for brittle fracture" [25], its capability to de- 
scribe a real physical system that fails under external load 
is questionable; All the data on fracture surfaces mor- 
phology were obtained by in-plane fracture experiments, 
mainly in mode I, whereas the random fuse model is a 
discrete version of out-of-plane fracture (mode III), which 
is unstable experimentally [26]. While in-plane fracture 
modes entail the solution of the fully tensorial Lame elas- 
ticity theory, out-of-plane fracture entails the solution of 
scalar elasticity theory, described by the Laplace equa- 
tion. There is no a-priori reason to believe that the two 
problems are in the same universality class. Moreover, 
the random fuse model misrepresents another important 
feature of fracture experiments; when a material is exter- 
nally loaded intrinsic defects can lead to the nucleation 
of a crack everywhere in the system. Since this process 
occurs in an uncontrolled way, it is an experimental con- 
vention to introduce a small crack in the sample prior 
to loading. Typically, the stress concentration near the 



tip of such a crack is so overwhelming that damage far 
from this region cannot develop significantly even in het- 
erogenous materials. This is in contrast to the random 
fuse model, excluding situations of very narrow disor- 
ders, where damage is accumulated throughout the sys- 
tem prior to the critical phase of macroscopic failure in 
which the final rupture line develops. In spite of the in- 
adequacy of the random fuse model to faithfully describe 
real physical systems, it might be useful to understand 
how positive correlations are built in the model, lead- 
ing to H > 1/2. Recently [25] it was proposed that a 
correlated gradient percolation process is responsible for 
the value of the roughness exponent measured in this 
model. Here we do not give the argument in full detail, 
but focus on the main conceptual aspect of this proposi- 
tion: correlations in the model are crucially dependent on 
the damage accumulated in the system before the critical 
phase of the formation of the final rupture line. In this 
view, the positive correlations observed in the final rup- 
ture line is a manifestation, via a coalescence process, of 
the long range positive correlations that are built in the 
system during the earlier stages of the dynamics. There- 
fore, this explanation rely heavily on the presence of a 
broad enough threshold distribution that is responsible 
for correlations in the precursory phase of the dynamics. 
On the other hand, other works [27, 28] support the view 
that correlations in the damage accumulated prior to the 
critical phase of macroscopic failure are negligible. If this 
view is valid then the origin of self-affine crack roughness 
in the random fuse model does not depend on whether 
there is "strong" or "weak" disorder since correlations 
are built in the system only at the final stage of macro- 
scopic failure. We think that this view represents better 
the typical experimental situation where stress concen- 
tration near the pre-existing crack dominates the failure 
process. This view will be further developed below, see 
Sec. IV. 

A second group of lattice models originated from the 
central force model [22] . In this model the lattice is con- 
nected by elastic springs that are freely rotating around 
the nodes. This model is much closer to real fracture 
since its continuum limit is the tensorial Lame elasticity. 
Unfortunately, no systematic studies of the roughness ex- 
ponent for this model were reported. A third group of lat- 
tice models is the random beam models. In these models 
the lattice is composed of elastic beams that are rigidly 
connected to the nodes. Here the bond bending elastic- 
ity is taken into account by considering the local rota- 
tional degrees of freedom at each node in addition to the 
translational ones that are described by the central force 
model. The roughness exponent in 1-1-1 dimensions for 
this model was found to be « 0.86 [29]. Therefore, this 
model seems to be unrelated to the experimental find- 
ings. The origin of this discrepancy is probably the fact 
that Lame elasticity provides an appropriate description 
of materials on a large enough length scale. To conclude, 
we propose that none of the lattice models found in the 
literature are capable of describing 1-|-1 dimensional in- 
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plane fracture and to faithfully represent the experiments 
in the field. 



B. Continuum Models 

The basic issue arising in continuum models of crack 
propagation is how to predict where the crack goes under 
the action of a given stress field. One can use symmetry 
principles to derive an equation of motion for the crack 
tip. Adapting Eq. (15) to two-dimensional quasi-static 
fracture, the stress tensor field attains the universal form 



Using these expressions Eq. (23) can be written as 



aij{r,e) = 



= S^/W, (21) 



where I and II denote the mode I (tensile) and mode II 

(shear) parts of the stress tensor field. As was explained 
above, the r~^/^ singularity is not physical; neverthe- 
less in many cases there exists a cut-off distance from 
the crack tip (a one-dimensional front), above which the 
stress field is dominated by the universal form. There- 
fore, the non-universal parts, the so-called stress intensity 
factors and K^^ are very important physical quanti- 
ties. 

In order to derive an equation of motion for the crack 
tip consider a curved crack in a two dimensional medium. 
The position of the crack tip is r*'^ in some convenient 
coordinate system. Let us denote the unit tangent and 
the unit normal to the curved crack at its tip by t and 
n respectively and assume that the crack tip fields are 
characterized by the stress intensity factors and ifjj . 
A discrete symmetry that the crack tip dynamics should 
remain invariant under its operation is n — > — n. Un- 
der this symmetry the relevant quantities in the problem 
transform as follows: (i) (ii) K-^-^ —Kji (iii) 

V V. Now one can write down the most general first 
order difi^erential equations that are invariant under this 
symmetry. The first equation is just the kinematic rela- 
tion for the rate of crack growth 



(22) 



The second one describes the rate of crack tip rotation 



^^=-f{K,,Kl,v)K,,h 



(23) 



where f{K^,Kf^,v) is a positive material function. These 

equations were originally derived in [30] and were shown 
to be consistent with various experimental observations 
in [31, 32]. To further simplify Eq. (23) we rewrite the 
local crack tip directions in terms of the angle 6 that the 
unit tangent makes with the x-axis as 



t = cos 6* X -I- sin y 
h = — sin ^ X -|- cos 6 y 



de 



-f{K,,Kl)K,, 



(25) 



In order to use this equation, consider a long 
crack propagating quasi-statically in a two dimensional 
medium under mode I conditions. If the medium were 
homogenous, we would have K^^{h(x)} = 0, and ac- 
cording to Eq. (25) the crack would extend along a 
straight path. Real material possess however some inho- 
mogeneities which induces fluctuations h{x) in the crack 
path. The global symmetry is broken locally, leading to 
K^^{h{x)} ^ 0. ifjj{/i(a;)} is a functional of the fluctua- 
tions h{x), typically including a long-range contribution, 
reflecting the long-range nature of the elastic interactions 
and a local contribution. If we restrict ourselves to small 
fluctuations h{x) we obtain 9 « d^h. Furthermore, by 
considering steady state propagation such that dt vdx 
and representing the effect of material disorder by a noise 
term '&{x), we arrive at the following Langevin type equa- 
tion 



7'^ 



,h{x) = -K^p{h{x)} - dix) 



(26) 



where the superscript (1) denotes the fact the K^^{h{x)} 
is considered to first order in the small fluctuation h{x) 
and the disorder is modeled as uncorrelated noise 



{^{x)^{x')) ~ S{x - x') 



(27) 



(24) 



The mode II stress intensity factor, K-^-^{h{x)}, was cal- 
culated in first order perturbation theory in [33], yielding 

2 V TT ^-x' 

(28) 

where the superscript (0) denotes quantity of the unper- 
turbed solution. By considering a configuration in which 
CTra (x') is practically constant, a scaling analysis of Eq. 
(26) results in a logarithmic roughness of the rupture 
line. A similar treatment as the one presented here was 
developed in [34] where higher order terms were shown 
to be irrelevant in the renormalization group sense. It is 
concluded that this approach does not produce any non- 
trivial roughness exponent. The reason for this failure is, 
in our opinion, the fact that the crack growth takes place 
at the tip which is controlled completely by the stress 
intensity factor K^^. We will show that there are obser- 
vations that motivate the view that the crack growth is a 
discrete process taking place via the nucleation of damage 
at a finite distance R ahead of the crack tip. This length 
R characterizes the scale of "spikiness" of the crack path; 
thus, when the growth step takes place at a distance R 
ahead of the crack tip the irregularities of the crack path 
( "spikiness" ) change the structure of the stress field and a 
different long range quantity controls the growth process. 
Moreover, the irregularities of the crack path are difficult 
to handle using perturbative approaches. The effect of 
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this finite length scale on the crack growth is further dis- 
cussed in the next section. We conclude this section by 
noting that the origin of the long range positive corre- 
lations in 1+1 dimensional fracture and its multiscaling 
nature are not yet explained by the existing models. 



IV. A MODEL 

In this section we review a recently introduced sta- 
tistical model of crack surfaces morphology in 1+1 di- 
mensional fracture which takes into account the insights 
gained from both experiments and previous theoretical 
approaches. It is common to divide the failure of mate- 
rials into two categories: ductile and brittle. In brittle 
fracture materials break when subjected to large enough 
stresses without large plastic deformations taking place. 
In ductile fracture, on the other hand, materials undergo 
significant plastic flow before experiencing any local fail- 
ure. In this case a crack evolves by the nucleation, growth 
and coalescence of voids. A recent experiment [35] mo- 
tivated the view that quasi-static brittle fracture occurs, 
similarly to ductile fracture, by the nucleation of damage 
ahead of the crack tip. It is not clear yet if this pro- 
cess necessarily involves plastic deformations. The length 
scale of damage nucleation R in traditional "brittle" and 
"ductile" materials is completely different. Nevertheless, 
we suggest that this behavior is generic irrespective of 
whether the physics behind this length scale is related to 
plastic deformation or to existing material disorder. 

Therefore, a statistical theory of crack surfaces mor- 
phology in 1+1 dimensions should include the following 
components: (i) the exact solution of the elasticity prob- 
lem in the presence of irregular crack geometries (ii) a 
growth law in which the evolution of the crack is con- 
trolled by the nucleation of voids at a finite distance R 
ahead of the crack tip (iii) an appropriate description of 
the effect of material disorder. 



A. The elasticity problem in the presence of 
irregular crack geometries 

In this section we review a non-perturbative approach 
to the calculation of the stress field in two space dimen- 
sions for an arbitrarily shaped crack, based on conformal 
mapping [36]. We start with the quasi-static version of 
Eq. (10), obtained by dropping the inertial term on the 
right hand side 



dxi 



= 0. 



(29) 



For in-plane modes of fractures in two dimensions one 
introduces the Airy stress potential U{x,y) such that 



dy 



1 — 



2 )"xi/ 



dxdy 



(30) 



Thus the set of Eq. (29), after simple manipulations, 
translate to a bi-Laplace equation for the Airy stress po- 
tential U{x, y) [3] 



V^V2C/(a;,y) = , 



(31) 



with the prescribed boundary conditions on the crack 
and on the external boundaries of the material. At this 
point we choose to focus on the case of uniform remote 
loadings and traction-free crack boundaries. This choice, 
although not the most general, is of great interest and 
will serve to elucidate our method. Other solutions may 
be obtained by superposition. Thus, the boundary con- 
ditions at infinity, for the two in-plane symmetry modes 
of fracture, are presented as 

o'xxioo) = ■,(7yy(co) = (Too -..(Jxyico) = Mode 1(32) 
o'xx(oo) = ; cTyy (oo) = ; axy (oo) = CToo Mode II . 

(33) 

In addition, the free boundary conditions on the crack 
are expressed as 



i(s) = cryn{s) = 



(34) 



where s is the arc-length parametrization of the crack 
boundary and the subscript n denotes the outward nor- 
mal direction fi. 

The solution of the bi-Laplace equation can be written 
in terms of two analytic functions (p{z) and r]{z) as 



U{x,y)^n[z(f{^.)+T](z)] 



(35) 



In terms of these two analytic functions, using Eq. (30), 
the stress components are given by 



ayy{x,y) = ^[2^'{z) + z^"{z) + v"{z)] 
<Jxx{x,y) = ^2ip'{z)-zip"{z)-r]"{z)] 
a^y{x,y) = ^z^"{z) + r]"{z)]. 



(36) 



In order to compute the full stress field one should first 

formulate the boundary conditions in terms of the ana- 
lytic functions Lp{z) and 77(2) and to remove the gauge 
freedom in Eq. (35). The boundary conditions Eq. (34), 
using Eq. (30), can be rewritten as [37] 



ds 



dU .dU_ 

dx dy 



. 



(37) 



Note that we do not have enough boundary conditions to 
determine U{x,y) uniquely. In fact we can allow in Eq. 
(35) arbitrary transformations of the form 



(p ^ ip + iCz + 7 

+ i' = 'n' 



(38) 



where C is a real constant and 7 and 7 are complex 

constants. This provides five degrees of freedom in the 
definition of the Airy potential. Two of these freedoms 
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are removed by choosing the gauge in Eq. (37) according 
to 



— i— - 
dx dy 



on the boundary 



(39) 



It is important to stress that whatever the choice of the 
five freedoms, the stress tensor is unaffected; see [37] for 
an exhaustive discussion of this point. Computing Eq. 
(39) in terms of Eq. (35) we arrive at the boundary 
condition 



Lp{z) + zif'{z) + %p{z) = 



(40) 



To proceed we represent '^{z) and '(/'(^) in Laurent ex- 
pansion form: 



ip{z) 
i;{z) 



ipiz + V^o + ip-i/z + 'iJj-2/z^ + 



(41) 



This form is in agreement with the boundary conditions 
at infinity that disallow higher order terms in z. One 
freedom is now used to choose tpi to be real and two 
more freedoms will allow us later on to fix tpo- Then, 
using the boundary conditions (32) and (33), we find 



^ ^ ^ ^ ^^^^ J 

4 2 



= 



^1 



Mode II 



(42) 



As said above, the direct determination of the stress 
tensor for a given arbitrary shaped crack is difficult. To 
overcome the difficulty we perform an intermediate step 
of determining the conformal map from the exterior of the 
unit circle to the exterior of our given crack. Before our 
work the best available approach for such a task was the 
Schwartz-Cristoffel transformation. We have presented 
an alternative new approach for finding the wanted con- 
formal transformation, given in terms of a functional it- 
eration of fundamental conformal maps. The use of it- 
erated conformal maps was pioneered by Hastings and 
Levitov [38]; it was subsequently turned into a powerful 
tool for the study of fractal and fracture growth pat- 
terns [39-45]. In the next subsection we describe how, 
given a crack shape, to construct a conformal map from 
the complex w-plane to the physical z-plane such that 
the conformal map z = $(a;) maps the exterior of the 
unit circle in the w-plane to the exterior of the crack 
in the physical z-plane, after n directed growth steps. 
We draw the reader's attention to the fact that this 
method is more general than its application in this pa- 
per [45], and in fact we offer it as a superior method 
to the Schwartz-Cristoffel transformation, with hitherto 
undetermined potential applications in a variety of two- 
dimensional contexts. 



1. The conformal mapping 

The essential building block in the present application, 
as in all the applications of the method of iterated con- 
formal maps is the fundamental map (j)\^0 that maps the 



z-plane 





co-plane 




FIG. 7: Example of how to construct the conformal mapping 
along a line, see text for details. 



exterior circle onto the unit circle with a semi-circular 
bump of linear size \/X which is centered at the point 



This map reads 



00, A (w) = Vw 



(1 + A) 
2w 



(l + w) 



2 1 - A 
wl + X 



1/2 



The inverse mapping ipg^Q ^ is of the form 



Xz - yi + xjz^ - 1) 

1- (1 + A)z2 ' 



(43) 
1/2 

(44) 
(45) 



By composing this map with itself n times with a judi- 
cious choice of series {8k}k=i ^-^d {Xk}n=i will con- 
struct $("^(a;) that will map the exterior of the circle 
to the exterior of an arbitrary simply connected shape. 
To understand how to choose the two series {Ok}^^i and 
{Xk}n=i consider Fig. 7 and define the inverse map uj = 
X*^"^ (z). Assume now that we already have (w) and 

therefore also its analytic inverse X^""^^-^) after n — 1 
growth steps and we want to perform the next iteration. 
To construct $'•"-'(0;) we advance our mapping in the di- 
rection of a point z in the z-plane by adding a bump in 
the direction oi w = x'"~^''(-?) f^ie w-plane. The map 
$(")(w) is obtained as follows: 

$(")(c^) = ((/.e„,A„(c^)) . 

The value of 0„ is determined by 

0„ = arg[x("-iHz)] 



(46) 



(47) 
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The magnitude of the bump A„ is determined by requir- 
ing fixed size bumps in the z-plane. This means that 



An — 



Ao 



$(n-i)'(eie„)|2 



(48) 



We note here that it is not necessary in principle to have 
fixed size bumps in the physical domain. In fact, adap- 
tive size blimps could lead to improvements in the pre- 
cision and performance of our scheme. We consider here 
the fixed size scheme for the sake of simplicity and we 
will show that the accuracy obtained is sufficient for our 
purposes. Iterating the scheme described above we end 
up with a conformal map that is written in terms of an 
iteration over the fundamental maps (43): 



(49) 



For the sake of newcomers to the art of iterated conformal 
maps we stress that this iterative structure is abnormal, 
in the sense that the order of iterates in inverted with 
respect to standard dynamical systems. On the other 
hand the inverse mapping follows a standard iterative 
scheme 



(50) 



The algorithm is then described as follows; first we 
divide the curve into segments separated by points {zi}. 
The spatial extent of each segment is taken to be approx- 
imately \f\), in order to match the size of the bumps in 
the z-plane. Without loss of generality we can take one 
of these points to be at the center of coordinates and to 
be our starting point. From the starting point we now 
advance along the shape by mapping the next point Zi 
on the curve according to the scheme described above. 



Using the Laurent form of the conformal map, Eq. (51), 
the linear term as w ^ oo is determined by Eqs. (52). 
We therefore write 

(^(w) = ifiFiuj + ip-i/u) + + ■ ■ ■ , 

i){uj) = Vi-Fiw-|-V!o + V!-i/w + Vi-2/w^ + ... ,(53) 

where we used the last two freedoms to choose (pQ = 
—Foifo such that (po = 0. The boundary condition (40) 
is now read for the unit circle in the w-plane. Denoting 
e = e*^ and 



fe=i 



fe=0 



we write 



«(e) + |^«'(e)+^;(e) = /(e) 

$ (e) 



(54) 



(55) 



The function /(e) is a known function that contains all 
the coefficients that were determined so far: 



/(e) = -t^iFie - =^(^iFi 

$ (e) e 



(56) 



To solve the problem we need to compute the coefficients 
(pn and V'n- To this aim we first write [44] 



1>'(e) 



= 5:6.e^ 



(57) 



The function /(e) has also an expansion of the form 

00 

/(e) = ^/,eV (58) 



2. Solution in terms of conformal mappings 

The conformal map $^"^(0;) is constructed in n itera- 
tive steps. For the discussion below we do not need the n 

superscript and will denote simply ^(lu) = <!>'•"•' (w). This 
map is univalent [39], having the Laurent expansion form 



$(a;) = Fiuj + Fo + F_i/(j + F_2/a;^ + ■ 



(51) 



Any position z in the physical domain is mapped by 
x{z) = $"^(z) onto a position w in the mathematical do- 
main. This transformation does not immediately provide 
the solution as the bi-Laplace equation, in contrast to the 
Laplace equation, is not conformally invariant. Never- 
theless, the conformal mapping method can be extended 
to non-Laplacian problems. Wc begin by writing our un- 
known functions ip{z) and ip{z) in terms of the conformal 
map 



(52) 



In the discussion below we assume that the coefficients 6j 
and fi arc known. In order to compute these coefficients 
we need to Fourier transform the function <I>(e)/$' (e). 
This is the most expensive step in our solution. One 
needs to carefully evaluate the Fourier integrals between 
the branch cuts. Using the last two equations together 
with Eqs. (54) and (55) we obtain 

00 

<P-m - Ykb-m-k-W-k = f-m , m = 1, 2 • • -(59) 
/c=l 

oo 

V'-m - hm-k-W-k = fm ■ "1 = 0, 1, 2 • • • (60) 

fe=l 

These sets of linear equations are well posed. The coef- 
ficients (p-m can be calculated from equation (59) alone, 
and then they can be used to determine the coefficients 
'i()-m- This is in fact a proof that Eq. (55) determines 
the functions u{e) and w(e) together. This fact had been 
proven with some generality in [37]. 
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The calculation of the Laurent expansion form of ip{uj) 
and ip{u)) provides the solution of the problem in the w- 
plane. Still, one should express the derivatives of ip{z) 
and 77(2) in terms of ip{Lu) and ifjiuj) and the inverse map 
X{z) to obtain the solution in the physical z-plane. This 
is straight forward and yields 

f'{z) = 'f'ixiz)] x'{z) 

cp"{z) = <p"[x{z)] [x'{z)]^ + 0'[x{z)]x"{z) 

V"{z) = ^'{z) = i>'[x{z)] X'{z). (61) 

Upon substituting these relations into Eq. (36) one can 

calculate the full stress field for an arbitrarily shaped 
crack. The expression of the stress field in terms of the 
inverse conformal mapping is known for quite a long time 
although it is very limited as the conformal mapping and 
its inverse is rarely at hand. The central step of progress 
here is the conjunction of the novel functional iterative 
scheme for obtaining the inverse conformal mapping with 
the known result that expresses the stress field in terms 
of this inverse mapping. This method was shown in [36] 
to reproduce results for known regular crack geometries. 



B. A Crack Growth Law 

After developing the mathematics needed for the calcu- 
lation of the linear elastic stress field for irregular crack 
geometries, we turn now to a description of the crack 
growth law and the effect of disorder. As was explained 
above, a crucial aspect of the crack growth process is the 
finite length scale R in which damage nucleates. A sim- 
ple model [46, 47] for R can be developed by assuming 
the near crack tip zone to be properly described by the 
Huber-von Mises plasticity theory [48]. This theory fo- 
cuses on the deviatoric stress = aij — |Tr(cr)(5y and 
on its invariants. The second invariant, J2 = ^SijSij, cor- 
responds to the distortional energy. The material yields 
as the distortional energy exceeds a material-dependent 
threshold cr^ • The fact that we treat this threshold as a 
constant, independent of the state of deformation and its 
history, implies that we specialize for "perfect" plasticity. 
In two-dimensions this yield condition reads [48] 



al — crifT2 -I- cr| 2 
7: = f,. 



Here 0-1,2 are the principal stresses given by 



± 



vv 



(62) 



(63) 



In the purely linear-elastic solution the crack tip region 
is where high stresses arc concentrated. Perfect plasticity 
implies on the one hand that the tip is blunted and on 
the other hand that inside the plastic zone the Huber-von 
Mises criterion (62) is satisfied. The outer boundary of 
the plastic zone will be called below the "yield curve". 



and in polar coordinates around the crack tip will be de- 
noted TZ{9). Below we will compute the outer stress field 
exactly as was explained in the previous section. Using 
this field we can find the yield curve TZ{9). Such a typical 
curves are shown in the insets of Figs. 8 and 10. 

Whatever is the actual shape of the blunted tip its 
boundary cannot support stress. Together with Eq. (62) 
this implies that on the crack interface 



0-1 



CT2 = 0. 



(64) 



The typical scale R follows from the physics of the nu- 
cleation process. It is physically plausible that void for- 
mation is more susceptible to hydrostatic tension than to 
distortional stresses. We assume that void nucleation oc- 
curs where the hydrostatic tension P, P = ^Trcr, exceeds 
some threshold value Pc- The calculated the hydrostatic 
tension P along the yield curve is found to be significantly 
higher than its value on the crack surface P = (cf. 
Eq. (64)). On the other hand the linear-elastic solution 
implies a monotonically decreasing P outside the yield 
curve. We thus expect P to attain its maximum value 
between the blunted crack tip and the yield curve. This 
conclusion is fully supported by finite elements method 
calculations, cf. [49]. As the detailed elastic- plastic crack 
tip fields are not computed here, we use the outer elastic 
solution on the yield curve to determine the void nucle- 
ation position. We expect that this approximation should 
not affect the roughening exponents on large scales. The 
distance from the crack tip where P attains its maximal 
value is of the order of the size of the plastic zone. 

In this model the nucleation occurs when P exceeds a 
threshold Pc- The void will thus appear at a typical dis- 
tance R. A very clear demonstration of the appearance 
of the void near the boundary of the plastic zone is seen 
in the molecular dynamics simulations of [50] . Note that 
R is not a newly found length scale; it is the well known 
scale of the plastic zone [51]. Its identification with a 
length scale that is related to the scaling behavior of frac- 
ture surfaces is however new. This stems from the propo- 
sition that positive correlations appear only between the 
positions of nucleated voids. Below R one enters the 
regime of plastic processes whose theory is far from be- 
ing settled. We should also comment that it is possible 
that positive correlations appear even below the scale of 
the plastic zone since experiments indicate that several 
voids nucleate within the plastic zone [52, 53]. Finally, we 
note that the assumption of perfect plasticity, i.e. that 
<Ty is independent of the state of deformation and its 
history, is not true for real materials; usually is not 
sharply defined; it can increase with plastic deformations 
[48]. This phenomenon, known as "work-hardening" or 
"strain-hardening" is not taken into account in this sim- 
ple model, with the hope to be irrelevant for the crack 
morphology on large scales. 

Naturally, the precise location of the nucleating void 
will experience a high degree of stochasticity due to 
material disorder. A basic assumption behind our 
treatment of this stochasticity is that the material 
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disorder occurs on scales smaller than i?, allowing us 
to introduce a probability distribution function that 
is a function of the hydrostatic tension on the yield 
curve. Since we do not know from first principles the 
probability distribution for void formation, we consider 
in our model below two possible distribution functions. 
In all cases nucleation cannot occur if P < Pc. For 
P > Pc the void occurs with probability 

P oc P-Pc , (65) 
P oc exp[a(P - Pc)] - 1 . (66) 

In the exponential case we considered two different val- 
ues of a. In Fig. 8 wc show three such pdf's as they 
appear for a perfectly straight crack. We note that these 
distributions are symmetric about the forward direction. 
Nevertheless they have sufficient width to allow devia- 
tions from forward growth. These deviations will be re- 
sponsible later for the roughening of the crack. For com- 
parison examine also the pdf's for a general crack which 
are shown in Fig. 10. There the symmetry is lost: corre- 
lation to previous steps create a preference for the upward 
direction. 

Each growth step in our model is composed of two 
events. Firstly the material yields near the crack tip, 
creating a plastic zone with a void growing somewhere 
at the zone boundary. Secondly the crack tip and the 
void coalesce. We note that there is a separation of time 
scales between these two events. The first is slow enough 
to be governed by a quasi-static stress field. The second 
event occurs on a shorter time scale. It is clear that 
we forsake in the model any detailed description of the 
geometry on scales smaller than R. Any relevant scaling 
exponent that will be found in this model will refer to 
roughening on length scales larger than R. Therefore, 
the physical process in which the crack coalesces with 
the multiple voids ahead of it is substituted by a single 
void coalescence with the crack. 



C. Results 

Calculating the stress field around the crack, according 
to the method of iterated conformal maps, we can readily 
find the yield curve and the physical region in its vicinity 
where a void can be nucleated. Choosing with any one of 
the probability distributions described above, we use this 
site as a pointer that directs the crack tip. We then use 
the method of iterated conformal mappings to make a 
growth step to coalesce the tip with the void. Naturally 
the step sizes are of the order of R. We reiterate that 
this model forsakes the details of the void structure and 
all the length scales below R. Since wc arc making linear 
steps below R, we anticipate having an artificial scaling 
exponent iJ = 1 for scales smaller than R. This is clearly 
acceptable as long as wc are mainly interested in the 
scaling properties on scales larger than R. 
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FIG. 8: Panel (i): the tip of a straight crack and the yield 
curve in front of it. 

Panel (ii): three probability distribution functions calculated 
for the configuration in (i). The abscissa is Q, the angle mea- 
sured from the crack tip as seen in panel (i). The ordinate 
is the normalized probability (per unit 6) to grow in the Q 
direction. The distributions arc symmetric and wide enough 
to allow deviations from the forward direction. For all the 
curves = 6. For curve (a) p(S) oc exp[(P — P^)] — 1 and 
^ = 8^ for curve (b) p(6l) oc exp[0.2(P - P^)] - 1 and Ji, = 6 
and for curve (c) p(^) oc P — Pc and = 6. 



In Fig. 9 we present two typical cracks that were grown 
using this method. Both cracks were initiated from a 
straight crack of length 10000, representing the exper- 
imental paradigm of introducing a notch before load- 
ing the sample. The upper crack was grown using the 
broader exponential pdf of Fig. 10 curve (b). The lower 
crack was grown with the narrower pdf of Fig. 10 curve 
(a). Clearly, the upper crack exhibits stronger height 
fluctuations, as can be expected from the wider pdf and 
the choice of parameters. For the lower crack forward 
growth is much more preferred. In the upper crack the 
positive correlations between successive void nucleation 
and coalescence events can be seen even with the naked 
eye. This is precisely the property that we were after. 
A neat way to see this tendency is in the pdf's as they 
are computed on the yields crack for a typical, rather 
than straight, crack. In Fig. 10 we show these pdf's for 
the crack whose yield curve is shown in the upper panel. 
We see that now the symmetry of the pdf's is lost, and 
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FIG. 9: Two typical cracks generated with our model. Note 
the huge difference in scales of the abscissa and ordinate (the 
cracks are globally flat) and that the lower crack had been 
translated by -300. The upper crack exhibits two decades of 
self-affinc scaling with a Hurst exponent 0.64. The lower crack 
has smaller standard deviation and therefore a shorter scaling 
range. Nevertheless it appears that in its shorter scaling range 
it exhibits an exponent that is very close to the upper crack. 



positive values of 9 are preferred. This is the source of 
positive correlations that eventually gives rise to a non- 
trivial roughening exponent. 

This is born out by the measurements of the scaling 
properties of the fracture lines morphology that we dis- 
cuss next. 

We first restrict ourselves to a monoscaling analysis. 
Duo to the significant computational cost of the iterated 
conformal maps technique the numerical investigation of 
the growth model had a limited number of realizations of 
a few hundreds growth steps. As a result of the relative 
paucity of data, the structure functions defined in Eq. 
(17) would not converge well enough to provide reliable 
exponents. Comparing the various available methods for 
estimating roughness exponents [54] we decided to select 
the max-min method, which seems to give reliable results 
in our case. Therefore, one defines S{i) according to 

S{t) = {Max{y{x)}^^-^^^, - Mm{y(i)}^<-<^+,)^ . 

(67) 

For self-afRne graphs the Hurst exponent H, as was ex- 
plained in the introduction, is obtained via the scaling 
relation 



(68) 



In Fig. 11 we present a typical log- log plot of S{t) vs. 
I, in this case for the two cracks in Fig. 9 with power- 
law fits oi H = 0.64 and H = 0.68 respectively. Indeed, 
as anticipated from the visual observation of Fig. 9 the 
exponent is higher than 0.5. 

It turned out that all the cracks grown by our algo- 
rithm gave rise to scaling plots in which a scaling range 
with H = 0.66 ± 0.03 is clearly seen. When the pdf al- 
lowed for a sizeable standard deviation, the cracks gave 
a very nice scaling plot with at least two decades of clear 
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FIG. 10: Panel (i): the tip of a "rough" crack and the yield 
curve in front of it. Panel (ii): three probability distribution 
functions calculated for the configuration in panel (i). The 
abscissa is ^, the angle measured from the crack tip as seen in 
panel (i). The ordinate is the normalized probability (per unit 
&) to grow in the B direction. The pdf 's are those used in Fig. 
8, using the same parameters. Note the upward preference in 
all the pdf 's due to the broken symmetry. 



anomalous scaling. When the standard deviation was 
small, the scaling range was more meager, as seen in 
Fig. 11. It is interesting to stress that the anomalous 
scaling exponent appears insensitive to the pdf used, al- 
though the extent of the scaling range clearly depended 
on the pdf. We note that our measured scaling exponents 
are very close to the exponents observed in other two- 
dimensional experiments [9-12]. In addition the value of 
R does not effect the scaling properties of a crack, i.e. it 
doesn't seem to matter how long the step is, so long as a 
wide distribution of angles is allowed. 

We now turn to a multiscaling analysis of the morphol- 
ogy of fracture lines. For this purpose, we define 

S^{i) ^ (|Maar{/i(i)},<,<,^,-Mm{/z(5-)},<,<,+, |") . 

(69) 

The scaling exponents Q"^' are defined in analogous way 
to Eq. (18) by 

S„(A^)~A«'"'5„(£) . (70) 
Note that Sx{i) is just S{i) defined in Eq. (67) and C^^^ 
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FIG. 11: Calculation of the roughening exponent H. The 
slopes of the dotted lines are 0.64 for the upper plot (curve 
a) and 0.68 for the lower (curve b). Note that the initial 
scaling with slope 1 is relevant for length scales smaller than 
R. This scaling is an artifact, resulting from our algorithm 
that connects the crack tip to a void by a straight line. 



is just H defined in Eq. (68). 




FIG. 12: The spectrum ("^"^ as a function of the moment order 
n for rupture lines in the model of [46, 47]. The function is 
fitted to the form C^"' = nH — n^X and the parameters H and 
A are given. The errors in the estimation of these parameters 
reflect both the variance between different realizations and fit 
quality. The linear plot n(^'^' is added for comparison. 

The resulting exponents C^"-' for the cracks generated 
by this model are shown in Fig. 12. For the low or- 
ders moments (here we are limited by the paucity of 
data to n < 5) one again fits a quadratic function, with 
H = 0.66 ± 0.03 and A = 0.023 ± 0.003. The errors in the 
estimation of these parameters reflect both the variance 
between different realizations and the fit quality. The 



n dependence of the exponents C^"-* and the values of 
the fitting parameters are in agreement with the exper- 
imental ones. Since there is nothing in the model that 
is specific for the physics of paper, it appears that mul- 
tiscaling is a generic property of the fracture process, at 
least in 1-1-1 dimensions. 

To sum up, we have developed a model that shows that 
long-range elastic interactions, the appearance of a finite 
length scale in the growth process and material disorder 
are able to reproduce the scaling properties associated 
with the morphology of 1-1-1 dimensional fracture. It 
should be stressed that a deeper understanding of the 
nature of the long-range correlations induced by geomet- 
rical irregularities of the crack is still lacking. It would 
be nice to specify which physical quantity controls the 
growth at a finite length scale R away from the crack 
tip. Moreover, it is important to study this model with 
different local crack tip physics to better understand the 
degree of universality and the role of the length scale R. 



V. SCALING PROPERTIES IN 2+1 
DIMENSIONAL FRACTURE 

In previous sections we have considered the scaling 
properties of l-fl dimensional rupture lines. In this sec- 
tion we extend the discussion to include fracture surfaces 
which are 2-1-1 dimensional graphs. In 2-1-1 dimensions 
one denotes the graph as h(r) and considers again the 
structure function S2{£), which is now a two dimensional 
function 



(71) 



where angular brackets denote an average over all r. 
Initially no attention was paid to the fact that the 
isotropy in the fracture plane is broken due to initial 
conditions that lead to a preferred propagation direc- 
tion and the statement of [7] was that the structure 
function is a homogeneous function of its arguments, 
S2{Xe.) ^ A«*''S'2(£), as implied by Eq. (16). In fact such 
a statement is tenable only if the fracture process and the 
material are isotropic. Usually the crack propagates pre- 
dominantly in one direction (say x) and the vector £ de- 
fines an angle 6 with respect to x, 9 = cos~^ {x ■ (.) . There 
is no reason why the scaling exponent (^^^\ if it exists at 
all, should not depend on the angle 6. Indeed, in the later 
work that followed [7] this problem was recognized and 
scaling exponents were sought for one dimensional cuts 
through S2{i), typically parallel and orthogonal to the 
direction of the crack propagation. Besides the obvious 
meaning of 'parallel' and 'orthogonal' to x, no reason was 
ever given why these particular directions are expected 
to provide clean scaling properties. We will argue below 
that in general such one dimensional cuts exhibit a mix- 
ture of scaling exponents with amplitudes that depend on 
the angle 9, where 6 = and 6 = tt/2 are not special. To 
better understand the way crack propagation directional- 
ity affects scaling isotropy, imagine a fracture experiment 
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FIG. 13: Sketch of a hypothetical fracture experiments ar- 
ranged to allow a crack to develop in an isotropic fashion, i.e. 
with all angles Q being statistically equivalent. On the left, the 
full three dimensional experiment is shown, where the tensile 
axis is along z and a circular cavity is in the xy plane. On 
the right, a magnified version of the circular cavity in the xy 
plane is shown. 



in which an initial circular cavity is made to propagate by 
a tensile load such that the crack edge remains circular 
on the average, without any preferred propagation direc- 
tion in the plane normal to the load, see Fig. 13. From 
the point of view of the scaling properties of the rough 
fracture surface that is left behind the advancing crack, 
such an experiment is the analogue of homogenous and 
isotropic turbulence in nonlinear fluid mechanics. Nor- 
mal experiments in both fracture and turbulence involve 
symmetry breaking; the boundary conditions introduce 
anisotropy, making the discussion of scaling properties 
non-trivial. In turbulence it was shown how to disen- 
tangle the anisotropic contributions from the isotropic 
one by projecting the measured correlation and structure 
functions on the irreducible representations of the SO (3) 
symmetry group [55]. The scaling phenomena seen in the 
isotropic sector of anisotropic experiments are identical 
to those expected in the hypothetical experiment of ho- 
mogenous and isotropic turbulence. In this section we 
describe how a similar concept is introduced to the field 
of fracture: we will show that decomposing the height- 
height structure functions of fracture surfaces into the 
irreducible representations of the SO (2) symmetry group 
results in a simplification and rationalization of the scal- 
ing properties that is not totally dissimilar to the one 
obtained in turbulence. The scaling properties of the 
isotropic sector should be observable in principle in an 
experiment like the one shown in Fig. 13, which con- 
trary to turbulence may be performed in reality. 

Given an experimental surface h,(r) we first compute 
the second order structure function Eq. (71). The vec- 
tor i is associated with a norm t and an angle Q. By 
construction, the second order structure function is sym- 
metric under — > -I- tt. Accordingly, decomposing the 
structure functions into the irreducible representations 
of the SO (2) symmetry group results in summations over 
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FIG. 14: The raw fracture surface of the aluminum alloy ob- 
tained in Ref. [56]. 



even indices only: 

oo 



,i2m9 



(72) 



m— — oo 



Such a decomposition is deemed useful when each of the 
scalar functions a2m is itself a homogeneous function of 
its argument, characterized by an m dependent exponent: 



|a2™(A^)|^A«-|a2™(^)| , 



(73) 



where | • | stands for the norm of a complex number. For 
an isotropic fracture in an isotropic medium we expect 
0'2m{^) = for all m 7^ 0. In usual mode I experiments 
in which the crack propagates along the x direction and 
the tensile load is in the normal direction, there should be 
the same physics along lines with angles 6 and —9. This 
invariance under 9—^—0 implies that the arguments of 
all i2m(^) 7^ should be or tt. In reality this invariance 
might not hold on large length scales due to the paucity 
of data or due to some symmetry breaking process, see 
below. 

Our first experimental example was obtained [56] from 
a compact tension specimen made of 7475 aluminum alloy 
first precracked in fatigue and then broken under tension 
in mode I. The raw fracture surface and the second order 
structure function computed from it are shown in Figs. 
14 and 15 respectively. 

One sees the anisotropy of 52 (i) with the naked eye. To 
quantitatively characterize this anisotropy, the structure 
bmction was decomposed as in Eq. (72). The log- log 
plots of ao(^), 2|a2(i!)| and 2|a4(£)| are exhibited in Fig. 
16. 

By performing linear fit of the relevant range in the 
log-log plots we find the following exponents 



/■(2) 
^0 - 



1.32±0.07 



/-(2) 



1.45±0.08 



/-(2) 



2.1±0.1. (74) 



The implication is that at smaller length-scales the 

(2) 

smaller exponent Q should be dominant and vice versa. 
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X mm 

FIG. 15: Contour plot of the second order structure function 
of the surface shown in Fig. 14. 




FIG. 16; Log-log plots of the amplitudes ao{£), 2\a2{£)\ and 
2\a4{(.)\ vs. £ for the aluminum alloy. 



Indeed, examining again the contour plot in Fig. 15 one 
observes that at small scales the contours tend to ellipses 
of smaller eccentricity, whereas at larger scale the con- 
tours are ellipses with increasing eccentricity. 

The crucial test of this approach is whether one can re- 
construct the structure function S2{t,0) using the func- 
tional form of the irreducible representation and a mini- 
mal number of parameters. Indeed, at smaller values of £ 
the first two irreducible representations suffice. Writing 

S2{1, 0) « 8.30 £1-32 ^ 3 ^1-45 005(26* + tt) , (75) 

we compare in Fig. 17 the experimental data to Eq. (75) 
for £ = 5, 15 and 25mm. The excellent fit is obvi- 
ous. In fact, with four parameters (two amplitudes and 
two exponents) we can represent the structure function 
to within 1% in norm as long as ^ < 30mm. For 
larger values of £ the agreement decreases, and we need 
to employ the next irreducible representation. Adding 
0.026 cos(46' + 7r), we find the fit shown in Fig. 18 for 
£ — 35mm. Beyond these values the power-laws fits lose 
their credibility for this experimental data set. 

A second experimental example was obtained from the 
dynamic fracture of artificial rocks produced from car- 
bonatic aggregates cemented by epoxy [57] . The samples 
are plates of size 400 x 400 x 9 mm and the fracture sur- 
face was measured using a scanning laser profilometer. 
The analysis of the experimental data follows verbatim 
the first example. The plots of ao(^) and 2|a2(£)| are 
shown in Fig. 19. 

Fitting the plots we find 

Co^'= 0.86 ± 0.05 , = 1.93 ± 0.05 , cf ^= 1-93 ± 0.1 . 

(76) 

Two comments are in order. First, one should notice the 
non-universality of the scaling exponents as compared 
with the previous example (74). Second, the present sur- 
face does not satisfy a 9 ^ —9 symmetry. As the ex- 
periments generating the surface were full dynamic, the 
fracture front undergoes the well-known microbranching 
instability, resulting in directed "branch lines" that break 
the 6 —9 symmetry [58]. Due to the lack of symmetry 
the amplitudes of the coefhcient a™ can take any phase, 
not constrained to or tt as required by the — > — 
symmetry. The lack of symmetry is clearly obvious in 
the reconstruction of the structure function from the ir- 
reducible representations. In Fig. 20 we compare the 
experimental values of S2{£,d), for £ = 2 mm, to the 
expansion 

S2{£,e) « 0.025 ^° ®^-f 0.0016 £i-^3cos(26' + 2.09) 

+ 5.4 X 10"'' £^-^^ cos{49 - 0.17) , (77) 

The fit is satisfactory and the asymmetry in 6 is obvious. 

Taking the present two examples as representative, it 
appears that the S0(2) decomposition extracts pure scal- 
ing behavior in each sector, but that the scaling expo- 
nents are not universal, at least in the two experiments 
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FIG. 17: The experimental S2{£,0) for the aluminum alloy 
(circles) and the representation Eq. (75) (line), for £ = 5, 15 
and 25mm. 



discussed here. A possible explanation for the discrep- 
ancy in the exponents is that the latter experiment was 
fully dynamic whereas the former was quasi-static. Con- 
sidering cuts in S2{i, 0) along the 9 — Q and 9 = it direc- 
tions, the present approach predicts a mixture of scaling 
exponents rather than pure power laws, potentially lead- 
ing to spurious exponents. 

Finally, one should point out that the SO (2) decompo- 
sition is not expected to yield satisfactory results when 
the material itself is strongly anisotropic. As an example 
we considered fracture surfaces in wood. This is clearly 
an anisotropic medium due to the fiber structure and in- 
deed we found that along and across the fiber directions 
the scaling behavior appears credible, whereas the S0(2) 
decomposition failed altogether to reveal clean scaling 
properties in any sector. 

In summary, we propose that materials which can be 
fractured in an isotropic fashion, i.e materials having an 
isotropic structure, often have anisotropic fracture sur- 
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FIG. 18: The experimental S2{(.,0) for the aluminum alloy 
(circles) and the SO (2) expansion up to the third even order 
irreducible representation (line), for I — 35mm. 
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FIG. 19: Log-log plots of the amplitudes ao{l) and 2|a2(^) 
for the artificial rock. 




FIG. 20: Comparison of the experimental S2{£,d) (circles) 
and the S0(2) expansion up to the third even irreducible rep- 
resentation (line) for the artificial rock with £ — 2mm. 
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faces only because of the breaking of isotropy by the ini- 
tial conditions. In such cases it appears useful to analyze 
the anisotropic contributions as "corrections to scaling" 
beyond the isotropic sector, which is always there, with 
a leading scaling exponent. The analysis reveals non- 
universality in the scaling exponents, a finding that calls 
for further future study and assessment, including the in- 
teresting question of the possible existence of universal- 
ity classes. On the practical side, we have demonstrated 
that the full information concerning the two dimensional 
structure function can be efficiently parameterized by a 
few amplitudes and scaling exponents. The reader should 
note that here we dealt only with second order structure 
functions. In analogy to turbulence it may be possible 
to decompose any higher order structure function into 
S0(2) irreducible representations [59]. This may reveal 
additional interesting scaling properties such as the phe- 
nomenon of multiscaling discussed in Sec. II. Finally, we 
would like to emphasize the great interest in the proposed 
isotropic fracture experiment and the measurement of the 
roughness exponent in such an experiment. If indeed this 
scaling exponent were identical to the exponent of the 
isotropic sector in a standard experiment, this would sig- 
nificantly strengthen the theoretical interest in the pro- 
posed approach. 

VI. SUMMARY 

The main thrust of this paper is that a careful study of 
the scaling properties of the graphs representing fracture 



surfaces can lead to a better understanding of the physics 
of fracture. We showed that in 1+1 dimensions fracture 
lines are multiscaling, a property that is not reproduced 
in a number of traditional models of fracture. Rather, 
one needs to consider possible failures of linear elasticity 
in the vicinity of the crack tip, to form a typical scale 
ahead of the crack tip characterizing the growth steps 
in the crack propagation. Also in 2+1 dimensions we 
find, after decomposing the statistical objects into their 
S0(2) irreducible representations, a host of anomalous 
exponents in the various sectors that cannot be under- 
stood without a careful rendering of the physics of frac- 
ture. We thus propose that future research should focus 
on the mechanisms for the failure of linear elasticity and 
how the physics discovered manifests itself in the scaling 
properties of the graphs of the fracture surfaces. 
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